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It is known that tlie digital waveguide (DW) method for solving the wave equation numerically 
on a grid can be manipulated into the form of the standard finite-difference time-domain (FDTD) 
method (also known as the "leapfrog" recursion). This paper derives a simple rule for going in the 
other direction, that is, converting the state variables of the FDTD recursion to corresponding wave 
variables in a DW simulation. Since boundary conditions and initial values are more intuitively 
transparent in the DW formulation, the simple means of converting back and forth can be useful in 
initializing and constructing boundaries for FDTD simulations. 
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I. INTRODUCTION 

The digital waveguide (DW) method has been used for 
many years to provide highly efficient algorithms for mu- 
sical sound synthesis based on physical models 0, ■ 
For a much longer time, finite-difference time-domain 
(FDTD) schemes have been used to simulate more gen- 
eral situations at generally higher cost |H 1^ • In re- 
cent years, there has been interest in relating these meth- 
ods to each other and in combining them for more gen- 
eral simulations. For example, modular hybrid methods 
have been devised which interconnect DW and FDTD 
simulations by means of a KW-pipe [T^ IT^. The basic 
idea of the KW-pipe adaptor is to convert the "Kirchoff 
variables" of the FDTD, such as string displacement, ve- 
locity, etc., to "wave variables" of the DW. The W vari- 
ables are regarded as the traveling-wave components of 
the K variables. 

In this paper, we present an alternative to the KW 
pipe. Instead of converting K variables to W variables, 
or vice versa, in the time domain, conversion formulas 
are derived with respect to the current state as a func- 
tion of spatial coordinates. As a result, it becomes sim- 
ple to convert any instantaneous state configuration from 
FDTD to DW form, or vice versa. Thus, instead of pro- 
viding the necessary time-domain filter to implement a 
KW pipe converting traveling-wave components to physi- 
cal displacement of a vibrating string, say, one may alter- 
natively set the displacement variables instantaneously to 
the values corresponding to a given set of traveling-wave 
components in the string model. Another benefit of the 
formulation is an exact physical interpretation of arbi- 
trary initial conditions and excitations in the K-variable 
FDTD method. Since the DW formulation is exact in 
principle (though bandlimited), while the FDTD is ap- 



proximate, even in principle, it can be argued that the 
true physical interpretation of the FDTD method is that 
given by the DW method. Since both methods gener- 
ate the same evolution of state from a common starting 
point, they may only differ in computational expense, nu- 
merical sensitivity, and in the details of supplying initial 
conditions and boundary conditions. 



II. IDEAL STRING WAVE EQUATION 

For definiteness, let's consider simulating the ideal vi- 
brating string, as shown in Fig. ^ 



string Tension 



e - Mass/Lengtii 



FIG. 1: The ideal vibrating string. 

The wave equation for the ideal (lossless, linear, flexi- 
ble) vibrating string depicted in Fig. ^ is given by 



Ky" = ey 



(1) 



where 



K = string tension 
e = linear mass density 
y = string displacement 
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y{t,x) 



y 

y = -§ty{t,x) 

y' = i^yit^^) 



and "=" means "is defined as." The wave equation is 
derived, e.g., in pol |. 



Finite Difference Time Domain Scheme 
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Using centered finite difference approximations (FDA) 
for the second-order partial derivatives, we obtain a finite 
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difference scheme for the ideal wave equation |l9ll3]|: 

.... ^ ^ y{t + T,x)~2yit,x) + yit-T,x) 
yit^x) w (2) 

2/ ~ (3) 

where T is the time samphng interval, and X is a spatial 
sampling interval. 

Substituting the FDA into the wave equation, choosing 
X = cT, where c = \jKlt is sound speed (normalized to 
c = 1 below), and sampling at times t = nT and x = mX , 
we obtain the following explicit finite diS'erence scheme 
for the string displacement: 

y{n + 1, m) = y{n, m + 1) + y{n, m — 1) — y{n — 1, m) (4) 

where the sampling intervals T and X have been normal- 
ized to 1. To initialize the recursion at time n = 0, past 
values are needed for all m (all points along the string) at 
time instants n = — 1 and n = —2. Then the string posi- 
tion may be computed for all m by ^ for n = 0,l,2,.... 
This has been called the FDTD or leapfrog finite differ- 
ence scheme 0- 



B. Digital Waveguide Scheme 

We now derive the digital waveguide formulation by 
sampling the traveling-wave solution to the wave equa- 
tion. It is easily checked that the lossless ID wave equa- 
tion Ky" ~ eij is solved by any string shape ywhich 
travels to the left or right with speed c = -y/ K/e ■ De- 
note right-going traveling waves in general by yr(t — x/c) 
and left-going traveling waves by j/; (i -f x/c) , where yr 
and yi are assumed twice-differentiable. Then, as is 
well known, the general class of solutions to the lossless, 
one-dimensional, second-order wave equation can be ex- 
pressed as 



yit,x) ^yr(t-^^ +yi + ^) ■ 
Sampling these traveling-wave solutions yields 



(5) 



y{nT, mX) = yr{nT — mX/c) + yi{nT + mX/c) 
= yr[{n - m)T] + yiiin + m)T] 

= y^ {n — m) + y^ {n + m) (6) 

where a superscript denotes a "right-going" 

traveling-wave component, and "— " denotes propagation 
to the "left" . This notation is similar to that used for 
acoustic-tube modeling of speech . 

Figure 12 shows a signal flow diagram for the computa- 
tional model of which is often called a digital wave- 
guide model (for the ideal string in this case) [2^ Is^l ■ 
Note that, by the sampling theorem, it is an exact model 
so long as the initial conditions and any ongoing additive 
excitations are bandlimited to less than half the temporal 
sampling rate fs = 1/T Appendix G]. 
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FIG. 2: Digital simulation of the ideal, lossless wave- 
guide with observation points at s = and x — 3X = 
3cT. (The symbol "z"!" denotes a one-sample delay.) 



Note also that the position along the string, Xm = 
mX = mcT meters, is laid out from left to right in the 
diagram, giving a physical interpretation to the horizon- 
tal direction in the diagram, even though spatial samples 
have been eliminated from explicit consideration. (The 
arguments of y"*" and y~ have physical units of time.) 

The left- and right-going traveling wave components 
are summed to produce a physical output according to 



y{nT, mX) = y^ {n ~ m) + y {n + m) 



(7) 



In Fig. 121 "transverse displacement outputs" have been 
arbitrarily placed at a; = and x — 3X. The diagram 
is similar to that of well known ladder and lattice digital 
filter structures ^3 1 except for the delays along the upper 
rail, the absence of scattering junctions, and the direct 
physical interpretation. 



C. FDTD and DW Equivalence 

The FDTD and DW recursions both compute time up- 
dates by forming fixed linear combinations of past state. 
As a result, each can be described in "state-space form" 
|29l Appendix E] by a constant matrix operator, the 
"state transition matrix" , which multiplies the state vec- 
tor at the current time to produce the state vector for the 
next time step. The FDTD operator propagates K vari- 
ables while the DW operator propagates W variables. We 
may show equivalence by (1) defining a one-to-one trans- 
formation which will convert K variables to W variables 
or vice versa, and (2) showing that given any common 
initial state for both schemes, the state transition matri- 
ces compute the same next state in both cases. 

The next section shows that the linear transformation 
from W to K variables. 



y(n, m) ^ y^ {n — m) + y {n + m), 



(8) 



for all n and m, sets up a one-to-one linear transfor- 
mation between the K and W variables. Assuming this 
holds, it only remains to be shown that the DW and 
FDTD schemes preserve mapping after a state tran- 
sition from one time to the next. While this has been 
shown previously |27j| . we repeat the derivation here for 
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completeness. We also provide a state-space analysis 
reaching the same conclusion in fJ3 

From Fig. El it is clear that the DW scheme preserves 
mapping © by definition. For the FDTD scheme, we 
expand the right-hand of 10} using (jSJ and verify that 
the left-hand side also satisfies the map, i.e., that y{n + 
1, m) ^ y^{n + 1 ~ m) + y~ {n + \ -\~ m) holds: 

y{n + 1, m) = m + 1) + m — 1) — y{n — 1, m) 
= — m — 1) + + m + 1) 

— m + 1) + + m — 1) 

—y^{n — m — 1) — + m — 1) 

= y^(n + m + 1) + i/+(n — m + 1) 
= 1) - ?7i] + 1) +m] 

= + 1, m) 

Since the DW method propagates sampled (bandlimited) 
solutions to the ideal wave equation without error, it fol- 
lows that the FDTD method docs the same thing, despite 
the relatively crude approximations made in (j2J). In par- 
ticular, it is known that the FDA introduces artificial 
damping when applied to first order partial derivatives 
arising in lumped, mass-spring systems [30l |. 

The equivalence of the DW and FDTD state transi- 
tions extends readily to the DW mesh [sol Is^ which is 
essentially a lattice-work of DWs for simulating mem- 
branes and volumes. The equivalence is more important 
in higher dimensions because the FDTD formulation re- 
quires less computations per node than the DW approach 
in higher dimensions (see for some quantitative com- 
parisons). 

Even in one dimension, the DW and finite-difference 
methods have unique advantages in particular situations 
[T^ . and as a result they are often combined together 
to form a hybrid traveling- wave/physical- variable simu- 
lation 0,|lll|ll|ll|llt2l|2|. 



III. STATE TRANSFORMATIONS 



model, all state variables arc defined as belonging to the 
same time n, as shown in Fig. |31 
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FIG. 3; DW flow diagram. 

From and referring to the notation defined in 
Fig. 121 we may write the conversion from W to K vari- 
ables as 

yn,m+l = Vn.rn+l + Vn.m+l 

+ 



yn.ni+l ^~ Vrt-in—l 



(9) 



where the last equality follows from the traveling-wave 
behavior (see Fig.O. 




In previous work, time-domain adaptors (digital fil- 
ters) converting between K variables and W variables 
have been devised 0|. In this section, an alternative 
approach is proposed. Mapping (jSJ gives us an imme- 
diate conversion from W to K state variables, so all we 
need now is the inverse map for any time n. This is com- 
plicated by the fact that non-local spatial dependencies 
can go indefinitely in one direction along the string, as we 
will see below. We will proceed by first writing down the 
conversion from W to K variables in matrix form, which 
is easy to do, and then invert that matrix. For simplicity, 
we will consider the case of an infinitely long string. 

To initialize a K variable simulation for starting at time 
n -|- 1, we need initial spatial samples at all positions m 
for two successive times n — 1 and n. From this state 
specification, the FDTD scheme Q can compute y{n + 
l,m) for all to, and so on for increasing n. In the DW 



FIG. 4: Stencil of the FDTD scheme. 

Figure |3| shows the so-called "stencil" of the FDTD 
scheme. The larger circles indicate the state at time n 
which can be used to compute the state at time n -\- 
1. The filled and unfilled circles indicate membership 
in one of two interleaved grids Q. To see why there 
arc two interleaved grids, note that when to is even, the 
update for yn+i.m depends only on odd to from time n 
and even to from time n—1. Since the two W components 
of yn-i,m. are converted to two W components at time 
n in we have that the update for yn+i,m depends 
only on W components from time n and positions to ± 1. 
Moving to the next position update, for y„+i^rn+i, the 
state used is independent of that used for yn+i.m, and the 
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W components used are from positions m and m-\-2. As the state- variable transformation separately for even and 
a result of these observations, we see that we may write odd m, e.g.^ 



yn,m — l 
l,m 

yn— l,m+4 
yn,7n+5 



1 1 

1 1 

1 1 

1 1 

1 

1 







1 

1 
1 



Vn^ra — l 
Vn^ra — \ 
Vn.m+l 
Vn.m+l 
Vn ,m+3 
,m+3 

Vn ,771 + 5 
Vn ,771 + 5 



(10) 



Denote the linear transformation operator by T and the 
K and W state vectors by and x^^ , respectively. Then 
(fTn|l can be restated as 

Tx^. (11) 

The operator T can be recognized as the Toeplitz op- 
erator associated with the linear, shift-invariant filter 
H{z) = 1 + z^^. While the present context is not a 



simple convolution since x^/ is not a simple time series, 
the inverse of T corresponds to the Toeplitz operator as- 
sociated with 

H{z) = — = l-z-^ + z-^-z-^ + -- - . 
1 + z ^ 

Therefore, we may easily write down the inverted trans- 
formation: 



y+ 



Vn , m — 1 
yn,m+l 
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yn,m+3 
Vn ,7n+5 

,771 + 5 
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±1 



Un — l.m 

y 71, 771 + 1 
^77-1,777+2 

^71,771 + 3 
^77-1, 771 + 4 

^77 , 771 + 5 



(12) 



The case of the finite string is identical to that of the 
infinite string when the matrix T is simply "cropped" 
to a finite square size (leaving an isolated 1 in the lower 
right corner); in such cases, as given above is simply 
cropped to the same size, retaining its upper triangular 
±1 structure. Another interesting set of cases is obtained 
by inserting a 1 in the lower-left corner of the cropped 
T matrix to make it circulant; in these cases, the M x 
M matrix contains ±1/2 in every position for even 
M, and is singular for odd M (when there is one zero 



eigenvalue) . 

IV. EXAMPLES 
A. Localized Displacement Excitations 

Whenever two adjacent components of are initial- 
ized with equal amplitude, only a single M^- variable will 
be affected. For example, the initial conditions (for time 
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1) 



?yn,m— 1 1 
yn-l,m = 1 

will initialize only y,7r?j.-i' ^ solitary left-going pulse of 
amplitude 1 at time n = 0, as can be seen from (|12ll by 
adding the leftmost columns explicitly written for T^^. 
Similarly, the initialization 



yn— l,?n — 2 1 
2/n,m — 1 1 

gives rise to an isolated right-going pulse y, 



+ 

71 , m - 



1' 



sponding to the leftmost column of plus the first 
column on the left not explicitly written in (|12|l . The su- 
perposition of these two examples corresponds to a phys- 
ical impulsive excitation at time and position to — 1: 



yn—l.m — 2 1 
2/n,?n— 1 2 
2/n — 1 , m — 1 



(13) 



Thus, the impulse starts out with amplitude 2 at time 
and position m— 1, and afterwards, impulses of amplitude 
1 propagate away to the left and right along the string. 

In summary, we see that to excite a single sample of 
displacement traveling in a single-direction, we must ex- 
cite equally a pair of adjacent colums in . This corre- 
sponds to equally weighted excitation of K- variable pairs 
the form (y„,m, 2/n-i,m±i)- 

Note that these examples involved only one of the two 
interleaved computational grids. Shifting over an odd 
number of spatial samples to the left or right would in- 
volve the other grid, as would shifting time forward or 
backward an odd number of samples. 



B. Localized Velocity Excitations 

Initial velocity excitations are straightforward in the 
DW paradigm, but can be less intuitive in the FDTD 
domain. It is well known that velocity in a displacement- 
wave DW simulation is determined by the difference of 
the right- and left-going waves [2^. Specifically, initial 
velocity waves can be computed from from initial dis- 
placement waves by spatially differentiating to ob- 
tain traveling slope waves y'^, multiplying by minus the 
tension K to obtain force waves, and finally dividing by 
the wave impedance R = \fK~t to obtain velocity waves: 

/■+ 

+ 1+ J 

w = -cy R 



f- 



(14) 



v~ {n+m) . (A more direct derivation can be based on dif- 
ferentiating |(SJ) with respect to x and solving for velocity 
traveling- wave components, considering left- and right- 
going cases separately at first, and arguing the general 
ease by superposition.) 

We can see from (|12|l that such asymmetry can be 
caused by unequal weighting of yn.m and yn,rn±i- For 
example, the initialization 

yn-l,m+l = +1 
Vn—l.Ta 1 

corresponds to an impulse velocity excitation at position 
TO -I- 1/2. In this case, both interleaved grids are excited. 



C. More General Velocity Excitations 

From H12(l . it is clear that initializing any single K vari- 
able yn.rn corrcsponds to the initialization of an infinite 
number of W variables y^,„ and y,7Tn- That is, a single 
K variable yn,m corrcsponds to only a single column of 
T^^ for only one of the interleaved grids. For example, 
referring to H12|) . initializing the K variable yn-i.m to -1 
at time n (with all other yn,m intialized to 0) corresponds 
to the W-variable initialization 

<m-(2,,+ l) = +1' M = 0,1,2,- •• 
Vm-(2M+1) = M = 0,1,2,- •• 

with all other W variables being initialized to zero. In 
view of earlier remarks, this corresponds to an impulsive 
velocity excitation on only one of the two subgrids. A 
schematic depiction from /i = TO— 5toTO-|-5of the W 
variables at time n is as follows: 
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(15) 



Below the solid line is the sum of the left- and right- 
going traveling- wave components, i.e., the corresponding 
K variables at time n. The vertical lines divide positions 
/i = TO — 1 and fj, = m. The initial displacement is zero 
everywhere at time n, consistent with an initial velocity 
excitation. At times ly = ?i + 1, n + 2, rt + 3, n + 4, the 
configuration evolves as follows: 
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where c = y/K/e denotes sound speed. The initial string 
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velocity at each point is then v{nT, mX) = v^{n — to) -|- 
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The sequence [...,1,0,1,0,1,...] consists of a dc (zero- 
frequency) component with amphtude 1/2, plus a sam- 
pled sinusoid of amplitude 1/2 oscillating at half the sam- 
pling rate fg = l/T. The dc component is physically 
correct for an initial velocity point-excitation (a spread- 
ing square pulse on the string is expected). However, 
the component at /s/2 is usually regarded as an artifact 
of the finite difference scheme. From the DW interpreta- 
tion of the FDTD scheme, which is an exact, bandlimitcd 
physical interpretation, we see that physically the com- 
ponent at /s/2 comes about from initializing velocity on 
only one of the two interleaved subgrids of the FDTD 
scheme. Any asymmetry in the excitation of the two in- 
terleaved grids will result in excitation of this frequency 
component. 

Due to the independent interleaved subgrids in the 
FDTD algorithm, it is nearly always non-physical to ex- 
cite only one of them, as the above example makes clear. 
It is analogous to illuminating only every other pixel in 
a digital image. However, joint excitation of both grids 
may be accomplished either by exciting adjacent spatial 
samples at the same time, or the same spatial sample at 
successive times instants. 

In addition to the W components being non-local, they 
can demand a larger dynamic range than the K variables. 
For example, if the entire semi-infinite string for m < 
is initialized with velocity 2X/T, the initial displacement 
traveling- wave components look as follows: 



(20) 
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and the variables evolve forward in time as follows: 



... 7 


6 


5 


4 


3 


2 


1 











••• 




••• -5 


-4 


-3 


-2 


-1 

















••• 


(21) 


••• 2 


2 


2 


2 


2 


2 


1 











••• 





... 8 
... „4 


7 

-3 


6 
-2 


5 4 3 
-10 


2 1 ••• 
••• 


... 4 


4 


4 


4 4 3 


2 1 ••• 


... 9 


8 


7 


6 5 4 


3 2 1 ••• 


... _3 


-2 


-1 





••• 


... 6 


6 


6 


6 5 4 


3 2 1 ••• 



(22) 



(23) 



Thus, the left semi-infinite string moves upward at a con- 
stant velocity of 2, while a ramp spreads out to the left 
and right of position fi = m at speed c, as expected phys- 
ically. By pO|l , the corresponding initial FDTD state for 
this case is 

yn,M = 0, ^j. eZ 

yn—l,ni—l 1 1 



Vn-l. 



/i < TO — 1, 



where Z denotes the set of all integers. While the FDTD 
excitation is also not local, of course, it is bounded for 
all fi. 

Since the traveling-wave components of initial veloc- 
ity excitations are generally non-local in a displacement- 
based simulation, as illustrated in the preceding exam- 
ples, it is often preferable to use velocity waves (or force 
waves) in the first place |33 • 

Another reason to prefer force or velocity waves is that 
displacement inputs are inherently impulsive. To see why 
this is so, consider that any physically correct driving in- 
put must effectively exert some finite force on the string, 
and this force is free to change arbitrarily over time. The 
"equivalent circuit" of the infinitely long string at the 
driving point is a "dashpot" having real, positive resis- 
tance 2R = The applied force f{t) can be divided 
by 2R to obtain the velocity v{t) of the string driving 
point, and this velocity is free to vary arbitrarily over 
time, proportional to the applied force. However, this ve- 
locity must be time-integrated to obtain a displacement 
y(t). Therefore, there can be no instantaneous displace- 
ment response to a finite driving force. In other words, 
any instantaneous effect of an input driving signal on 
an output displacement sample is non-physical except 
in the case of a massless system. Infinite force is re- 
quired to move the string instantaneously. In sampled 
displacement simulations, we must interpret displace- 
ment changes as resulting from time-integration over a 
sampling period. As the sampling rate increases, any 
physically meaningful displacement driving signal must 
converge to zero. 



D. Additive Inputs 

Instead of initial conditions, ongoing input signals can 
be defined analogously. For example, feeding an input 
signal Un into the FDTD via 



Un.m—l yn,m — 1 ^ ^n — 1 

yn,ra — yn,m ^~ 2u^ 
yn,m-(-l yn,m+l ^"^n— 1 



(24) 



corresponds to physically driving a single sample of string 
displacement at position to. This is the spatially dis- 
tributed alternative to the temporally distributed solu- 
tion of feeding an input to a single displacement sample 
via the filter H{z) = 1 — z"^ as discussed in [T^ . 
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E. Physical Interpretation of H{z) — 1 — z ^ 

As shown above, driving a single displacement sample 
Tjn.m in the FDTD corresponds to driving a velocity in- 
put at position m on two alternating subgrids over time. 
Therefore, the filter H{z) = 1 — acts as the filter 
H{z) = 1 — on either subgrid alone — a first-order 
difference. Since displacement is being simulated, ve- 
locity inputs must be numerically integrated. The first- 
order difference can be seen as canceling this integration, 
thereby converting a velocity input to a displacement in- 
put, as in (Pl|l . 



V. STATE SPACE FORMULATION 

In this section, we will summarize and extend the 
above discussion by means of a siate space analysis |l2j |. 

A. FDTD State Space Model 

Let x_i^{n) denote the FDTD state for one of the two 
subgrids at time n, as defined by The other subgrid 
is handled identically and will not be considered explic- 
itly. In fact, the other subgrid can be dropped altogether 
to obtain a half-rate, staggered grid scheme |3lllll|. How- 
ever, boundary conditions and input signals will couple 
the subgrids, in general. To land on the same subgrid 
after a state update, it is necessary to advance time by 
two samples instead of one. The state-space model for 
one subgrid of the FDTD model of the ideal string may 
then be written as 

Xf^{n + 2) = Ak Xj^{n) -^Bk u{n + 2) 

y{n) = CKXKin). (25) 

To avoid the issue of boundary conditions for now, we 
will continue working with the infinitely long string. As 
a result, the state vector Xj^{n) denotes a point in a space 
of countably infinite dimensionality. A proper treatment 
of this case would be in terms of operator theory [2ll| . 
However, matrix notation is also clear and will be used 
below. Boundary conditions are taken up in W CI 

When there is a general input signal vector u(n), it is 
necessary to augment the input matrix B/f to accomo- 
date contributions over both time steps. This is because 
inputs to positions m ± 1 at time n + l affect position m 
at time n + 2. Henceforth, we assume B/f and u have 
been augmented in this way. Thus, if there are q input 
signals v{n) = [vi{n)], i ~ 1, . . . ,q, driving the full string 
state through weights (3^ = [/3m. i], m € Z, the vector 
u(n) = is of dimension 2q x 1: 



u{n + 2) = 



vin 



1 and u is 2 X 1. The matrix B/^ weights these inputs 
before they are added to the state vector for time n -I- 2, 
and its entries are derived in terms of the (3m, i coefficients 
below. 

Ck forms the output signal as an arbitrary linear com- 
bination of states. To obtain the usual displacement out- 
put for the subgrid, Ck is the matrix formed from the 
identity matrix by deleting every other row, thereby re- 
taining all displacement samples at time n and discarding 
all displacement samples at time n — 1 in the state vector 
x^(n): 



771 — 2 




yn,m 




yn,m+2 




yn,Tn+4 




y{n) 



1 

1 

1 

1 



yn,m-2 
yn—l,m — l 

yn,m 
J/n-l^m+l 

yn,m+2 
yn— l,m+3 

yn,rn+4: 



The state transition matrix A/^ may be obtained by first 
performing a one-step time update. 



yn+2,'. 



^71+1, 771— 1 yn+1, 771+1 yn,m 

+ 13' v{n + 2), 



and then expanding the two n -I- 1 terms in terms of the 
state at time n: 

2/ri+l,m-l = yn.,,n-2 + yn.m " yri-l,m-l + 0^^_^ll{n + 1) 
yn+l,m+l = yn,m + 2/n,m+2 - yn-l,m+l + P^^^^Hin + 1) 

(26) 

The intra-grid state update for even m is then given by 

yn+2,m 

= yn,m-2 " yn-l,m-l + ^n,™ ^ 2/ri-l/iTi+l + yn,m+2 

+ (3^v{n + 2) + {P ,+f3^yv{n+l) 

—m — ^ ^ ^—771- 1 !—m-\-l-' — ^ ' 

[1,-1,1,-1,1] 

yn — 1,771 — 1 

y 71,771 

yn-l,m-|-l 

y?!, 771 + 2 

,yn — 1,771+3 , 



!—m ^—771 — 1 ■ — 771+ 1 ^ 



v{n + 2) 
v{n + 1) 



(27) 



When there is only one physical input, as is typically 
assumed for plucked, struck, and bowed strings, then q ~ 



For odd m, the update in (|26|l is used. Thus, every 
other row of A^, for time n -f 2, consists of the vector 
[1,-1, 1, —1, 1] preceded and followed by zeros. Succes- 
sive rows for time n + 2 are shifted right two places. The 
rows for time n-\-l consist of the vector [1,-1,1] aligned 
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similarly: 



yn + l,m-l 




yn + 2,m 















1 


-1 


1 














1 


-1 


1 


-1 


1 














1 


-1 


1 














1 


-1 


1 


-1 


1 



yn,m-2 
y 71 — 1 . m — 1 
yn.m 

yn~l.m + l 

yn,m+2 

2/n-l,m+3 



From (|27|l we also see that the input matrix is given 
as defined in the following expression: 



2/ri+l,Tn-l 

2/n+2,m 

yn+l,m+l 

J/n+2,m+2 

2/n+l,m+3 

2/n+2,m+4 

yn+l,m+5 





■ — m 



T 

■m+2 



tLm+4 





■ — m+o ^ — m+5 
i-m+5 



i;(n + 2) 
v(ri + 1) 

u(n+2) 



B. DW State Space Model 



As discussed in ijllll the traveling-wave decomposition 
(O defines a linear transformation from the DW 
state to the FDTD state: 



(28) 



Since T is invertible, it qualifies as a linear transfor- 
mation for performing a change of coordinates for the 
state space. Substituting l|28(l into the FDTD state space 
model (|25() gives 

Tx^y{n + 2) = AKTxyir{n)+BKu{n + 2) (29) 
y{n) = CifTx^(n). (30) 

Multiplying through (^5)) by gives a new state-space 
representation of the same dynamic system which we will 
show is in fact the DW model of Fig. 



Xy^^{n + 2) 
y{n) 



■ Bwu{n + 2) 



(31) 



where 



Aw = T-^Aa-T 



A 



To verify that the DW model derived in this manner is 
the computation diagrammed in Fig. |3 we may write 
down the state transition matrix for one subgrid from 
the figure to obtain the permutation matrix Ayy, 



Vn+2,m-2 




yn+2,m-2 




Vn+2,m 




Vn+2,m 




Vn+2,m+2 




Vn+2,m+2 









1000000000 
0000010000 
0010000000 
0000000100 
0000100000 
0000000001 



and displacement output matrix C 



Vn^m — i 
yn,m — i 
Vn,m-2 
?/n,m-2 
yn,m 

yn,m + 2 
yn,m + 2 
yn,m+4 
yn,m-\-4 



(33) 



w- 







Un^ra 








yn,'m-\-4: 











1 1 

1 1 

1 1 

1 1 



yn,m— 2 
2/n,m-2 

Vn^ra 
2/n,m+2 
Vn,m+2 
Vn,m+4 



1. DW Displacement Inputs 
We define general DW inputs as follows: 

vXm = Vn-l.m-l + {1^^ ^{n) 

-1 + ilZfuin) 



Vn-l. 



(34) 
(35) 



The mth 2q x 2 block of the input matrix Biy driv- 
ing state components {yn+2 rrn Vn+2 mV ^'^'^ multiplying 
[v_{n + 2)'^,v{n + 1)'^]'^ is then given by 



(7+)^ (7+ )^ 

ii'V (7" V 



(36) 



(32) 



Typically, input signals are injected equally to the left 
and right along the string, in which case 

7^ = 7 = 7 . 

— m —rn —rn 

Physically, this corresponds to applied forces at a single, 
non-moving, string position over time. The state update 
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with this simplification appears as 



FDTD input-weights according to 



yn+2,' 
J/n+2,' 



2w("+2) 



T 

y 

-m 



7 7 
_Lm — m — 1 



T 

y 



7 7 



u(n + 2) 
v{n + 1) 

i(n+2) 



2/n+l,m-l 



Vn + 1 
yn+2 



rn+1 
m+2 



Tm ^" T?n— 1 T?n— 1 ^ Tm— 1 



Tm ^ 7t7 



7m-l + 7m+l 



Tm ^ Tm+1 T?n+lT7nH-l 
7m+2 + 7m+2 7m+l + 7m+3 



27(71 + 2) 
W(71 + 1) 



Note that if there are no inputs driving the adjacent sub- 
grid (7 , 

grid scheme, the input reduces to 



7 , = 0), such as in a half-rate staggered 



The left- and right-going input-weight superscripts in- 
dicate the origin of each coefficient. Setting 7+ = 
results in 



x^^{n + 2) = Awx, 



7 

— m 

,T 

-m+2 
,T 

-m+2 



v{n + 2). 



To show that the directly obtained FDTD and DW 
state-space models correspond to the same dynamic sys- 
tem, it remains to verify that Ayy = T^^Aj^T. It is 
somewhat easier to show that 



TA 



w 



AkT 



1000010000 
0010010000 
0010000100 
0000100100 



A straightforward calculation verifies that the above 
identity holds, as expected. One can similarly verify 
Cw = Ck T, as expected. The relation B^^ = T^^ B/^ 
provides a recipe for translating any choice of input sig- 
nals for the FDTD model to equivalent inputs for the 
DW model, or vice versa. For example, in the scalar 
input case {q — 1), the DW input- weights 'Byy become 



B 



K 



7m+7m-l 27„_i 

27m 7m-l + 7'm+l 
7m + 7m+l 27,„+i 

27m+2 7m+l + 7m+3 



(37) 



Finally, when 7,„ = 1 and 7^ = for all fi =/= m, 
obtain the result familiar from (|24|l : 



we 



B 



K 



1 

2 
1 



Similarly, setting 7^ = for all fj, m + 1, the weight- 
ing pattern (1, 2, 1) appears in the second column, shifted 
down one row. Thus, Ba' in general (for physically sta- 
tionary displacement inputs) can be seen as the super- 
position of weight patterns (1, 2, 1) in the left column for 
even m, and the right column for odd m (the other sub- 
grid), where the 2 is aligned with the driven sample. This 
is the general collection of displacement inputs. 



2. DW Non-Displacement Inputs 

Since a displacement input at position m corre- 
sponds to symmetrically exciting the right- and left-going 
traveling- wave components and y", it is of interest 
to understand what it means to excite these components 
antisymmetrically. As discussed in j|IV CI an antisym- 
metric excitation of traveling-wave components can be 
interpreted as a velocity excitation. It was noted that 
localized velocity excitations in the FDTD generally cor- 
respond to non-localized velocity excitations in the DW, 
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and that velocity in the DW is proportional to the spa- 
tial derivative of the difference between the left-going and 
right-going traveling displacement-wave components (see 
(|14|l ). More generally, the antisymmetric component of 
displacement-wave excitation can be expressed in terms 
of any wave variable which is linearly independent rela- 
tive to displacement, such as acceleration, slope, force, 
momentum, and so on. Since the state space of a vibrat- 
ing string (and other mechanical systems) is traditionally 
taken to be position and velocity, it is perhaps most nat- 
ural to relate the antisymmetric excitation component to 
velocity. 

In practice, the simplest way to handle a velocity input 
Vmin) in a DW simulation is to first pass it through a 
first-order integrator of the form 



Hiz) 



1 



1 



1 



(38) 



to convert it to a displacement input. By the equivalence 
of the DW and FDTD models, this works equally well 
for the FDTD model. However, in view of ^IV CI this 
approach does not take full advantage of the ability of 
the FDTD scheme to provide localized velocity inputs for 
applications such as simulating a piano hammer strike. 
The FDTD provides such velocity inputs for "free" while 
the DW requires the external integrator H38|) . 

Note, by the way, that these "integrals" (both that 
done internally by the FDTD and that done by H38I) ') are 
merely sums over discrete time — not true integrals. As 
a result, they are exact only at dc (and also trivially at 
/s/2, where the output amplitude is zero). Discrete sums 
can also be considered exact integrators for impulse-train 
inputs — a point of view sometimes useful when interpret- 
ing simulation results. For normal bandlimited signals, 
discrete sums most accurately approximate integrals in a 
neighborhood of dc. The KW-pipe filter H{z) = 1 — z~'^ 
has analogous properties. 



3. Input Locality 

The DW state-space model is given in terms of the 
FDTD state-space model by (|32ll . The similarity trans- 
formation matrix T is bidiagonal, so that Ck and Cw = 
Cx T are both approximately diagonal when the out- 
put is string displacement for all m. However, since 
given in H12() is upper triangular, the input matrix 
= T~^Bx can replace sparse input matrices Ba' 
with only half-sparse "Bw, unless successive columns of 
are equally weighted, as discussed in illVI We can 
say that local K-variable excitations may correspond to 
non-local W- variable excitations. From H3()|) and 13 7|) . we 
see that displacement inputs are always local in both sys- 
tems. Therefore, local FDTD and non-local DW excita- 
tions can only occur when a variable dual to displacement 
is being excited, such as velocity. If the external integra- 
tor H38|l is used, all inputs arc ultimately displacement 
inputs, and the distinction disappears. 



C. Boundary Conditions 

The relations of the previous section do not hold ex- 
actly when the string length is finite. A finite-length 
string forces consideration of boundary conditions. In 
this section, we will introduce boundary conditions as 
perturbations of the state transition matrix. In addition, 
we will use the DW-FDTD equivalence to obtain phys- 
ically well behaved boundary conditions for the FDTD 
method. 

Consider an ideal vibrating string with AI = 8 spa- 
tial samples. This is a sufficiently large number to make 
clear most of the repeating patterns in the general case. 
Introducing boundary conditions is most straightforward 
in the DW paradigm. We therefore begin with the or- 
der 8 DW model, for which the state vector (for the 0th 
subgrid) will be 



Vnfi 
Vnfi 
2/n,2 



yn,2 

Vnfi 
Vnfi 



The displacement output matrix is given by 



1 1 

1 1 

1 1 

1 1 



and the input matrix B^k is an arbitrary M x 2q ma- 
trix. We will choose a scalar input signal u{n) driving 
the displacement of the second spatial sample with unit 
gain: 







1/2 1/2 

1/2 1/2 











B 



w 



The state transition matrix Aw is obtained by reduc- 
ing H33|l to finite order in some way, thereby introducing 
boundary conditions. 



1. Resistive Terminations 

Let's begin with simple "resistive" terminations at the 
string cndpoints, resulting in the reflection coefficient g 
at each end of the string, where |g| < 1 corresponds to 
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nonnegative (passive) termination resistances [30|- In- 
spection of Fig. O makes it clear that terminating the left 
endpoint may be accomplished by setting 



and the right termination corresponds to 

By allowing an additional two samples of round-trip delay 
in each endpoint reflectance (one sample in the chosen 
subgrid), we can implement these reflections within the 
state-transition matrix: 



"0 


m 

















0" 











1 














1 






































1 














1 






































1 














1 





























9r 






(39) 



The simplest choice of state transformation matrix T is 
obtained by cropping it to size M x M: 



1 1 
1 1 
1 



1 
1 



1 1 
1 
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1 





1 1 
1 



An advantage of this choice is that its inverse T"-'^ is 
similarly a simple cropping of the M = oo case. However, 
the corresponding FDTD system is not so elegant: 



A 



TAvi/T 



"0 


91 


-91 


hi 


~hi 


hi 


-hi 


hi - 


1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 

















1 


-1 


hr 


— hr 




















9r 


-9r_ 



where hi = \ + gi and hr = 1 + Qr- We see that 
the left FDTD termination is non-local for g ^ —1, 
while the right termination is local (to two adjacent spa- 
tial samples) for all g. This can be viewed as a con- 
sequence of having ordered the FDTD state variables 

[yn,m; ^/n — l,m+l 1 ■ ■ ■ 

] instead of [y„_i,„i 
Choosing the other ordering interchanges the endpoint 
behavior. Call these orderings Type I and Type II, re- 
spectively. Then T// ~ Tj ; that is, the similarity trans- 
formation matrix T is transposed when converting from 



Type I to Type II or vice versa. By anechoically coupling 
a Type I FDTD simulation on the right with a Type II 
simulation on the left, general resistive terminations may 
be obtained on both ends which are localized to two spa- 
tial samples. 

In nearly all musical sound synthesis applications, at 
least one of the string endpoints is modeled as rigidly 
clamped at the "nut". Therefore, since the FDTD, as de- 
fined here, most naturally provides a clamped endpoint 
on the left, with more general localized terminations pos- 
sible on the right, we will proceed with this case for sim- 
plicity in what follows. Thus, we set gi = —1 and obtain 



Aa' = 



"0 


-1 


1 

















1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 

















1 


-1 


1 + gr 


-1 - 9r 




















9r 


-9r 



2. Boundary Conditions as Perturbations 

To study the effect of boundary conditions on the state 
transition matrices Aw and A^-, it is convenient to write 
the terminated transition matrix as the sum of of the 

h 

"left-clamped" case Aw (for which gi ~ —1) plus a scries 
of one or more rank-one perturbations. For example, 
introducing a right termination with reflectance gr can 
be written 



Aw — Aw 



i-W 



^1,2 + 9rSs,7, 



(40) 



where 5ij is the M x M matrix containing a 1 in its 
(i, j)th entry, and zero elsewhere. (Following established 
convention, rows and columns in matrices are numbered 
from 1.) 

h 

In general, when i + j is odd, adding dij to Aw cor- 
responds to a connection from left-going waves to right- 
going waves, or vice versa (see Fig.O. When i is odd and 
j is even, the connection flows from the right-going to the 
left-going signal path, thus providing a termination (or 
partial termination) on the right. Left terminations flow 
from the bottom to the top rail in Fig. |21 and in such 
connections i is even and j is odd. The spatial sample 
numbers involved in the connection arc 2[(i — 1)/2J and 
2[(j — 1)/2J, where [xj denotes the greatest integer less 
than or equal to x. 

The rank-one perturbation of the DW transition ma- 
trix H40|) corresponds to the following rank-one perturba- 

h 

tion of the FDTD transition matrix Ak'- 



A 



Ak = Ak + 9^8,7 



12 



where 



^8,7 



0- 



0000000 
1 -1 
1 -1. 



.(41) 



In general, we have 



(42) 



Thus, the general rule is that 6ij transforms to a matrix 
A^j which is zero in all but two rows (or all but one 
row when i = 1). The nonzero rows are numbered i and 
i — 1 (or just i when i = 1), and they are identical, being 
zero in columns 1 : j — 1, and containing [1,-1,1,-1,...] 
starting in column j. 



3. Reactive Terminations 

In typical string models for virtual musical instru- 
ments, the "nut end" of the string is rigidly clamped 
while the "bridge end" is terminated in a passive re- 
flectance S{z). The condition for passivity of the re- 
flectance is simply that its gain be bounded by 1 at all 
frequencies |3fll |: 



S'(e^"^)|<l, VtjTe [-7r,7r). 



(43) 



A very simple case, used, for example, in the Karplus- 
Strong plucked-string algorithm, is the two-point-average 
filter: 



S{z) 



l + z- 



To impose this lowpass-filtered reflectance on the right 
in the chosen subgrid, we may form 

1 . 1 . 

Aw — Aw — 2 ^'^ ~ 2 

which results in the FDTD transition matrix 



which is still passive, since it obeys H43|) . but it does not 
have the desired amplitude response: Instead, it has a 
notch (gain of 0) at one-fourth the sampling rate, and 
the gain comes back up to 1 at half the sampling rate. 
In a full-rate scheme, the two-point-average fllter must 
straddle both subgrids. 

Another often-used string termination filter in digital 
waveguide models is specified by 30] 



s{n) = 
S{e^^^) = 



h I h 
4' 2' 4 

-7wT l + ft.cos(a;T) 
3 9 ■ 



where g G (0, 1) is an overall gain factor that affects the 
decay rate of all frequencies equally, while h £ (0, 1) con- 
trols the relative decay rate of low-frequencies and high 
frequencies. An advantage of this termination filter is 
that the delay is always one sample, for all frequencies 
and for all parameter settings; as a result, the tuning of 
the string is invariant with respect to termination filter- 
ing. In this case, the perturbation is 

Aw - Aw— Y'^(^^"5'^^)"f^(^^~3'^^)~X'^^^'^~"^'*^^ 

and, using 1)42(1 . the order A/ = 8 FDTD state transition 
matrix is given by 



Ak = 



where 






-1 


1 

















1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 











,91 


-91 


1+.92 


-1-52 


1+.93 


-1 - 33 








.91 


-91 


92 


-.92 


93 


-93 








A 


gh 














91 = 
















A 
















92 = 




91 












A 


gh 














.93 = 




f .92- 







Ak 



"0 


-1 


1 

















1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 

















1 


-1 


1 


-1 


1 

















1 


-1 


1 

















1/2 


-1/2 




















-1/2 


1/2 


-1 


-1 



This gives the desired filter in a half-rate, staggered grid 
case. In the full-rate case, the termination filter is really 



Siz) 



1 + z- 



The filtered termination examples of this section gen- 
eralize immediately to arbitrary finite-impulse response 
(FIR) termination filters S{z). Denote the impulse re- 
sponse of the termination filter by 

s{n) = [so, si, S2, . . . , sn], 

where the length N of the filter does not exceed Af/2. 
Due to the DW-FDTD equivalence, the general stability 
condition is stated very simply as 



\S{e^^^) 



N-l 

n=0 



< 1, VcjT € [-7r,7r). 
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^. Interior Scattering Junctions 

A so-called Kelly-Lochhaum scattering junction |l8ll3Cl| | 
can be introduced into the string at the fourth sample by 
the following perturbation 



Ak = AA' + (l-fc()A5,3 + fcrA5,8 + fc/A6,3 + (l-fcr)A 



6,8- 



Here, fc; denotes the reflection coefBcient "seen" from left 
to right, and fc,- is the reflectance of the junction from the 
right. When the scattering junction is caused by a change 
in string density or tension, we have = —ki. When it is 
caused by an externally imposed termination (such as a 
plectrum or piano-hammer touching the string), we have 
kj. ^ ki, and the reflectances may become filters instead 
of real values in [—1,1]. Energy conservation demands 
that the transmission coefficients be amplitude comple- 
mentary with respect to the reflection coefficients [33 ■ 

A single time-varying scattering junction provides a 
reasonable model for plucking, striking, or bowing a 
string at a point. Several adjacent scattering junctions 
can model a distributed interaction, such as a piano ham- 
mer, finger, or finite- width bow spanning several string 
samples. 

Note that scattering junctions separated by o ne sp atial 
sample (as typical in "digital wavcg uide filters" |30|) wifl 
couple the formerly independent subgrids. If scattering 
junctions are confined to one subgrid, they are separated 
by two samples of delay instead of one, resulting in round- 
trip transfer functions of the form H{z^) (as occurs in 
the digital waveguide mesh). In the context of a half- 
rate staggered- grid scheme, they can provide general IIR 
filtering in the form of a ladder digital filter p^l30l|. 



D. Lossy Vibration 

The DW and FDTD state-space models are equivalent 
with respect to lossy traveling- wave simulation. Figure [3 
shows the flow diagram for the case of simple attenuation 
by g per sample of wave propagation, where g e (0, 1] for 
a passive string. 

2/n,m— 1 yn,m ?/n,m+l 



9 

U0 



9 

-yn,n 

9 



yn,m—l yn,m Vn^m+l 

FIG. 5: DW flow diagram in the lossy case. 



The DW state update can be written in this case as 

Xy^[n + 2) = g'^Awx^,{n) -I- Bvyu(n -I- 2). 

where the loss associated with two time steps has been in- 
corporated into the chosen subgrid for physical accuracy. 
(The neglected subgrid may now be considered lossless.) 
In changing coordinates to the FDTD scheme, the gain 
factor g^ can remain factored out, yielding 



(n + 2) = g'^AKXjiin) + BKu{n + 2). 



When the input is zero after a particular time, such as 
in a plucked or struck string simulation, the losses can 
be implemented at the flnal output, and only when an 
output is required, e.g., 

y{n) = g"-yo{n) 

where yo(n) denotes the corresponding lossless simula- 
tion. When there is a general input signal, the state 
vector needs to be properly attenuated by losses. In the 
DW case, the losses can be lumped at two points per 
spatial input and output . 



E. State Space Summary 

We have seen that the DW and FDTD schemes cor- 
respond to state-space models which are related to each 
other by a simple change of coordinates (similarity trans- 
formation). It is well known that such systems exhibit 
the same transfer functions, have the same modes, and 
so on. In short, they are the same linear dynamic sys- 
tem. Differences may exist with respect to spatial locality 
of input signals, initial conditions, and boundary condi- 
tions. 

State-space analysis was used to translate initial con- 
ditions and boundary conditions from one case to the 
other. Passive terminations in the DW paradigm were 
translated to passive terminations for the FDTD scheme, 
and FDTD excitations were translated to the DW case 
in order to interpret them physically. 



VI. COMPUTATIONAL COMPLEXITY 

The DW model is more efficient in one dimension be- 
cause it can make use of delay lines to obtain an 0(1) 
computation per time sample |2^, whereas the FDTD 
scheme is 0{M) per sample (M being the number of 
spatial samples along the string) . There is apparently no 
known way to achieve 0(1) complexity for the FDTD 
scheme. In higher dimensions, i.e., when simulating 
membranes and volumes, the delay-line advantage dis- 
appears, and the FDTD scheme has the lower operation 
count (and memory storage requirements). 
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VII. SUMMARY 

An explicit linear transformation was derived for 
converting state variables of the finite-difference time- 
domain (FDTD) scheme to those of the digital waveguide 
(DW) scheme. The equivalence of the FDTD and DW 
state transitions was reviewed, and the proof of state- 
space equivalence was completed. Since the DW scheme 
is exact within its bandwidth (being a sampled traveling- 
wave scheme instead of a finite difference scheme), it 
can be put forth as the proper physical interpretation of 
the FDTD scheme, and consequently be used to provide 
physically accurate initial conditions and excitations for 
the FDTD method. For its part, the FDTD method pro- 
vides lower cost relative to the DW method in dimensions 
higher than one (for simulating membranes, volumes, and 
so on), and can be preferred in highly distributed non- 



linear string simulation applications. 

VIII. FUTURE WORK 

The simple state translation formulas derived here for 
the one-dimensional case do not extend simply to higher 
dimensions. While straightforward extensions to higher 
dimensions are presumed to exist, a simple and intu- 
itive result such as found here for the ID case could be 
more useful for initializing and driving FDTD mesh sim- 
ulations from a physical point of view. In particular, 
spatially localized initial conditions and boundary con- 
ditions in the DW framework should map to localized 
counterparts in the FDTD scheme. A generalization of 
the Toeplitz operator T having a known closed-form in- 
verse could be useful in higher dimensions. 
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